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Abstract 
We observe for the first time the process ce — nhe with data collected by the BESIII experi- 
ment. Significant signals are observed at the center-of-mass energy ys = 4.226 GeV, and the Born 
cross section is measured to be 53; + 2.7) pb. Evidence for nh; is observed at \/s = 4.358 GeV 
with a Born cross section of (10.0434 + 2.6) pb, and upper limits on the production cross section 


at other center-of-mass energies between 4.085 and 4.600 GeV are determined. 


PACS numbers: 13.25.Gv, 13.66.Bc, 14.40.Pq, 14.40.Rt 


201711.02267v1 


chinaXiv 


I. INTRODUCTION 


The spectroscopy of charmonium states below the open charm threshold is well es- 
tablished, but the situation above the threshold is more complicated. From the inclu- 
sive hadronic cross section in e*e^ annihilation, some vector charmonium states, (3770), 
(4040), (4160), (4415) are known with properties as expected in the quark model "i 
However, besides these states, several new vector states, namely the Y(4260), Y(4360) and 
Y (4660), have been discovered experimentally HI In addition, some new states with other 
quantum number configurations are also found in experiment, such as the X (3872), Z.(3900) 
and Z.(4020) states i lg. The common properties of these states are their relatively 
narrow width for decaying into a pair of charmed mesons, and their strong coupling to hid- 
den charm final states. Therefore, it is hard to explain all these resonances as charmonia 
and they are named ‘charmonium-like states’ collectively. Several unconventional explana- 
tions, such as hybrid charmonium 19] tetraquark bd p. hadronic molecule 25), 
diquarks ES SP kinematical effects have been suggested. See also Ref. 32, 
and references therein for a recent review. 

To understand the nature of these charmonium-like states, it is mandatory to investigate 
both open and hidden charm decays. Most of the observed vector charmonium-like states 
transit to spin-triplet charmonium states with large rate since the spin alignment of the c and 
c-quarks does not need to be changed between initial and final states. However, the spin-flip 
process ce — mh. has also been observed by the CLEO [34] and BESIII experiments 
e [35], and the large cross section exceeds theoretical expectations [sg]. Furthermore, two 
new structures have been reported ine*e^ > c *- he [33]. This may suggest the existence of 
hybrid charmonium states with a pair of cc in spin-singlet configuration which easily couples 
to an he final state. Consequently, searching for the process e*e^ — nh, will provide more 
information about the spin-flip transition, and the structures observed in ce — sh, 
may be observed also in the 7h. process. In addition, the transition Y(4S) — nh, has 
been observed in the bottomonium system [37]. The analogous process in the charmonium 
system is worth searching for to understand the dynamics in the 7 transition between heavy 
quarkonia. 

The CLEO Collaboration observed evidence of about 30 for ee — nh, based on 
586 pb! data taken at va = 4.17 GeV GA and the measured cross section is (4.7 42.2) pb. 
In comparison, BESIII has collected data samples of about 4.7 fb! in total at vs > 4.0 GeV. 
In this paper, a search is performed for the process e*e^ — nh, with he — yn, based on data 
samples collected with the BESIII detector at center-of-mass (c.m.) energies from 4.085 to 
4.600 GeV, as listed in Table [I| The integrated luminosities of these data samples are mea- 
sured by analyzing large-angle Bhabha scattering events with an uncertainty of 1.096 [38], 
and the c.m. energies are measured using the di-muon process [39]. In the analysis, ne is 
reconstructed with 16 hadronic final states: pp, 2(1*» ),2(K^ K-), K* K mm, ppn^m , 
Sra) OK 90074) KA. pon, eka, ROR wm omms RG, 
2(rta-)n, win nen, and 2(1*-- eraf, in which K$ is reconstructed from its m+ decay, 
and 7° and 7 from their yy final state. 


II. DETECTOR AND DATA SAMPLES 


BEPCII is a two-ring ete~ collider designed for a peak luminosity of 10? cm ?s^! at 


a beam current of 0.93 A per beam. The cylindrical core of the BESIII detector con- 
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sists of a helium-gas-based main drift chamber (MDC) for charged-particle tracking and 
particle identification (PID) through the specific energy loss dE’/dz, a plastic scintillator 
time-of-flight (TOF) system for additional PID, and a 6240-crystal CsI( TI) electromagnetic 
calorimeter (EMC) for electron identification and photon detection. These components 
are all enclosed in a superconducting solenoidal magnet providing a 1-T magnetic field. 
The solenoid is supported by an octagonal flux-return yoke instrumented with resistive- 
plate-counter muon detector modules interleaved with steel. The geometrical acceptance for 
charged tracks and photons is 93% of 47, and the resolutions for charged-track momentum 
at 1 GeV is 0.5%. The resolutions of photon energy in barrel and end-cap regions are 2.5% 
and 596, respectively. More details on the features and capabilities of BESIII are provided 
in Ref. : 

A Monte Carlo (MC) simulation is used to determine the detection efficiency and to esti- 
mate physics background. The detector response is modelled with a GEANT4-based "mi 
detector simulation package. Signal and background processes are generated with specialized 
models that have been packaged and customized for BESIII. 40,000 MC events are generated 
for each decay mode of 7, at each c.m. energy with KKMC and BESEVTGEN lu. The 
events are generated with an h, mass of 3525.28 MeV/c? and a width of 1.0MeV. The E1 
transition he — Yne is generated with an angular distribution of 1 + cos? 6*, where Hr is 
the angle of the E1 photon with respect to the he helicity direction in the he rest frame. 
Multi-body ne decays are generated uniformly in phase space. In order to study poten- 
tial backgrounds, inclusive MC samples with the same size as the data are produced at 
Vs = 4.23, 4.26 and 4.36 GeV. They are generated using KKMC, which includes the decay 
of Y (4260), ISR production of the vector charmonium states, charmed meson production, 
QED events, and continuum processes. The known decay modes of the resonances are gen- 
erated with BESEVTGEN with branching fractions set to the world_average values [45]. The 
remaining charmonium decays are generated with LUNDCHARM 4g], while other hadronic 
events are generated with PYTHIA 


III. EVENT SELECTION AND STUDY OF BACKGROUND 


According to the MC simulation of ee — nh, with he > yn. at V/s = 4.226 GeV, 
the energy of the photon emitted in the EI transition he — Yne is expected to be in the 
range (400, 600) MeV in the laboratory frame. Therefore, the signal event should have one 
E1 photon candidate with energy located in the expected region and one 7 candidate with 
recoil mass in the region of (3480,3600) MeV/c?. We define the 7 recoil mass Miei (17) 
as Miecoi(7])?c* = (Eem — Ey)? — [Pem — Pn c, where (Eem, Pem) and (Ep, Py) are the four- 
momenta of the e*e^ system and 7 in the ce rest frame. Since the E1 photon energy 
distribution in the laboratory frame will broaden with increasing c.m. energy, the energy 
window requirement is enlarged to (350, 650) MeV for the data sets collected at ys > 
4.416 GeV. The 7, candidate is reconstructed by the hadronic systems determined by the 
corresponding decay mode. The invariant mass of the hadronic systems is required to be 
within the mass range of (2940, 3020) MeV/c?. For the selected candidates, we apply a fit 
to the distribution of the 7 recoil mass to obtain the signal yield. 

Charged tracks in BESIII are reconstructed from MDC hits within a fiducial range of 
| cos 0| < 0.93, where 0 is the polar angle of the track. We require that the point of closest 
approach (POCA) to the interaction point (IP) is within 10cm in the beam direction and 
within 1cm in the plane perpendicular to the beam direction. A vertex fit constrains the 
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production vertex, which is determined run by run, and all the charged tracks to a common 
vertex. Since the K$ has a relatively long lifetime, it will travel a certain distance in the 
detector to the point where it decays into daughter particles. The requirements on the track 
POCA and the vertex fit mentioned above are therefore not applied to its daughter particles. 
The TOF and dE/dx information are combined to form PID confidence levels (C.L.) for 
the pion, kaon, and proton hypotheses; both PID and kinematic fit information is used to 
determine the particle type of each charged track, as discussed below. 

Electromagnetic showers are reconstructed by clustering EMC crystal energies. Efficiency 
and energy resolution are improved by including energy deposits in nearby TOF counters. 
A photon candidate is defined by showers detected with the EMC exceeding a threshold 
of 25 MeV in the barrel region (| cos0| < 0.8) or of 50 MeV in the end-cap region (0.86 < 
|cos0| < 0.92). Showers in the transition region between the barrel and the end-cap are 
excluded because of the poor reconstruction. Moreover, EMC cluster timing requirements 
are used to suppress electronic noise and energy deposits unrelated to the event. 

Candidates for 7° (7) mesons are reconstructed from pairs of photons with an invariant 
mass M (vy) satisfying |M(y7) — Mroml| < 15 MeV/c. A one-constraint (1C) kinematic 
fit with the M(^y) constrained to the 7° (y) nominal mass mpo (Mp) is performed 
to improve the energy resolution. We reconstruct K} — mtm candidates with pairs of 
oppositely charged tracks with an invariant mass in the mass range of |M(a7) — mx,| < 
20 MeV/c?. Here, mx, denotes the nominal mass of K9 [45]. A vertex fit constrains the 
charged tracks to a common decay vertex, and the corrected track parameters are used to 
calculate the invariant mass. To reject random zz" combinations, a kinematic constraint 
between the production and decay vertices, called a secondary-vertex fit, is employed [48], 
and the decay length is required to be more than twice the vertex resolution. 

The 7, candidate is reconstructed in its decay to one of the 16 decay modes mentioned 
earlier. After the above selection, a four-constraint (4C) kinematic fit is performed for each 
event imposing overall energy-momentum conservation, and the x26 is required to be less 
than 25 to suppress background events with different final states. If multiple 7. candidates 
are found in an event, only the one with the smallest x? = X46--X1o-- Xia - Xvertex 1$ retained, 
where vie is the x? of the 1C fit for 7° (7), x2, is the sum over all charged tracks of the 
x? of the PID hypotheses, and x2,4,, is the X? of the KÌ secondary-vertex fit. If more than 
one 7 candidate with recoil mass in the he signal region (3480 < Myecoi(7) < 3600 MeV/c?) 
is found, the one which leads to a mass of the 7, candidate closest to the ne nominal mass 
ma, is selected to reconstruct the ne- 

The requirement on x2c and mass (energy) windows for 7, ne and E1 photon reconstruc- 
tion are determined by maximizing the figure-of-merit, FOM = Ns/4/Ns + Ng, where Ng 
represents the number of signal events determined by MC simulation, and Nr represents 
the number of background events obtained from A, sidebands in the data sample. The 
cross section of e*e^ — nh, measured by CLEO m and the 7, branching ratios given by 
the Particle Data Group (PDG) are used to scale the number of signal events in the 
optimization. 

After applying all the criteria to the data sample taken at vs =4.226 GeV, the events 
cluster in the signal region in the two-dimensional distribution as shown in Fig. Ha). If the 
two-dimensional histogram is projected to each axis, clear rj. and he signals can be found in 
the expected regions as shown in Fig.[I(b) and (c). Meanwhile, no structure is observed in 
the events from the 7. (he) sideband regions. To further understand the background shape, 
events located in the 7 sideband regions are also investigated, which are shown by the green 
shaded area in Fig. [| (d) and are well described by a smooth distribution. 


T 


In addition, inclusive MC samples generated at va = 4.23 GeV are analyzed to study the 
background components. Here, the ratios among different components are fixed according 
to theoretical calculation or experimental measurements, except for the Bhabha process. 
A sample of 1.0 x 10’ Bhabha events (about 296 of the Bhabha events in real data) is 
generated with the BABAYAGA generator for background estimation. From this study, 
the dominant background sources are found to be continuum processes according to the 
MC truth information, while Y (4260) decays only give a small contribution to the total 
background. Most background events from resonance decays are z7J/w, wXco and open 
charm production. A similar conclusion can be drawn for data samples taken at other c.m. 
energies. From the study above, we conclude that the background shape in the 7 recoil mass 
can be described by a linear function. 
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FIG. 1: Mass spectrum obtained at \/s = 4.226 GeV. (a) The two-dimensional distribution of 
the invariant mass of the hadronic system and the recoil mass of 7; (b) mass of hadrons in he 
signal ([3.51, 3.55] GeV/c?) and sideband regions ([3.48, 3.50] GeV/c? and [3.56, 3.58] GeV /c?); 
(c) 7 recoil mass in Ne signal ([2.94, 3.02] GeV/c?) and sideband region ([2.87, 2.91] GeV/c? and 
[3.05,3.09] GeV /c?), and (d) yy recoil mass in 7) signal ([0.531,0.563] GeV/c?) and sideband regions 
([0.505, 0.521] GeV/c? and [0.573, 0.589] GeV/c?). For(b), (c), and (d), the dots with error bars rep- 
resent the distributions in the signal regions and the shaded histograms represent the distributions 
in the sidebands. 


IV. FIT TO THE RECOIL MASS OF 5 


To obtain the he yield for each ne decay channel, the 16 7 recoil mass distributions are 
fitted simultaneously using an unbinned maximum likelihood method. In the fit, the signal 
shape is determined by the MC simulation and the background shape is described by a linear 
function. The total signal yield of 16 channels is set to be Nobs, which is the common variable 
for all sub-samples and required to be positive. Nops x f; is the signal yield of the i-th channel. 
Here, f; refers to the weight factor f; = B;c;/ X €;B;, in which the B; denotes the branching 
fraction of 5, decays to the i-th final state and e; represents the corresponding efficiency. 
The efficiency for two-body ne decays is about 20%, for three- or four-body decays is about 
10% and for six-body decays it is about 6%. The signal and the background normalization 
for each mode are free parameters in the fit. The mode-by-mode and summed fit results 
are shown in Figs. P] and DI respectively. The x? per degree of freedom (dof) for this fit is 
x?^/dof = 17.2/15 = 1.15, where sparsely populated bins are combined so that there are at 
least 7 counts per bin in the y? calculation. The total signal yield is 41 +9 with a statistical 
significance of 5.8 o. 
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FIG. 2: Simultaneously fitted 7 recoil mass spectra in e* e^ — nhe, he > ^ne, Ne — Xi for the 16 
final states X; at ys —4.226 GeV. The dots with error bars represent the 7 recoil mass spectrum in 
data. 'The solid lines show the total fit function and the dashed lines are the background component 
of the fit. 


With the same method, evidence for e*e^ — nh, is found in the data sample taken at 
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FIG. 3: Sum of the simultaneous fits to 7 recoil mass spectra for all 16 ne decay modes at \/s = 
4.226 GeV. The dots with error bars represent the 7 recoil mass spectrum in data. The solid line 
shows the total fit function and the dashed line is the background component of the fit. The shaded 
histogram shows the events from the ne sidebands. 


ys = 4.358 GeV, as shown in Fig. Hl but no obvious signals are observed for the data sets 
taken at other c.m. energies. 


V. BORN CROSS SECTION MEASUREMENT 


The Born cross section is calculated using the following formula: 


Nobs 
lar uc RES CC C CEU IW RUE GE AR UE UE IP A A 1 
LOF ALA EBG > y) BO > 1) EB 0) 
Here, £ is the integrated luminosity of the data sample taken at each c.m. energy. (1+8) 
is the radiative correction factor, which is defined as 


f e(s(1 — z))F(z, s)dx 

c(s) : 
where F(z,s) is the radiator function, which is known from a QED calculation with an 
accuracy of 0.196 50]. Here, s is squared c.m. energy, and s(1 — x) is the squared c.m. 


aoP™ (ete > nhe) = 
(1 +ô) = 


(2) 
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FIG. 4: Sum of the simultaneous fits to 7 recoil mass spectra for all 16 ne decay modes at \/s = 
4.358 GeV. The dots with error bars represent the 7 recoil mass spectrum in data. The solid line 
shows the total fit function and the dashed line is the background component of the fit. The shaded 
histogram shows the events from 7j. sidebands. 


energy after emission of the ISR photons. o(s) is the energy dependent Born cross section 
in the range of [4.07, 4.6] GeV. Actually, the radiative correction depends on the Born cross 
section from the production threshold to the e*e^ collision energy, which is also what we 
want to measure in this analysis. Therefore, the final Born cross section is obtained in 
an iterative way. The efficiencies from a set of signal MC samples without any radiative 
correction are used to calculate a first approximation to the observed cross section. Then, by 
taking the observed cross sections as inputs, new MC samples are generated with radiative 
correction and the efficiencies as well as (1 + ô) are updated. After that, the cross sections 
can also be recalculated accordingly. The iterations are performed in this way until a stable 
result is obtained. The values of (1 + A) from the last iteration are shown in Table [I] 

The term |l + I|? is the vacuum-polarization (VP) correction factor, which includes 
leptonic and hadronic contributions. This factor is calculated with the package provided 
in Ref. [51]. The package provides leptonic and hadronic VP both in the space-like and 
time-like regions. For the leptonic VP the complete one- and two-loop results and the 
known high-energy approximation for the three-loop corrections are included. The hadronic 
contributions are given in tabulated form in the subroutine HARD5N [3]. The |1-- HJ? values 
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are also shown in Table 

Table[Jand Fig. [5]show the energy dependent Born cross sections from this measurement. 
Taking into account the CLEO measurement at ys = 4.17 GeV [34], the cross section 
from 4.085 ~ 4.600 GeV is parameterized as the coherent sum of three Breit-Wigner (BW) 
functions, as shown by the solid line in Fig. In the fit, the parameters of the BW 
around 4.36 GeV are fixed to those of the Y (4360) H while the parameters of the other 
two BW functions are left free in the fit. The fitted parameters of the free BW are: M, = 
(4204 3-6) MeV/c?, T, = (32422) MeV and M; = (4496 +26) MeV/c?, To = (1044-69) MeV, 
where the uncertainties are statistical. 


60 
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FIG. 5: Fit to the cross section of e* e^ — nh, as a function of c.m. energies. The square with 
error bar shows the measurement from CLEO , the dots with error bars refer to the results of 
this measurement, and the solid line shows the fit result with 3 coherent BW functions. 


VI. SYSTEMATIC UNCERTAINTIES 


In this section, the study of the systematic uncertainty for the cross section measurement 
at \/s = 4.226 GeV is described. The same method is applied to the other c.m. energies. 

The main contributions to the systematic uncertainties are from the luminosity mea- 
surement, the fit method, B(h. > yn.)B(n — yy), ISR correction, VP correction and 
X Bin, + Xi). The systematic uncertainties from different sources are listed in Table 
All sources are treated as uncorrelated, so the total systematic uncertainty is obtained by 
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TABLE I: Data sets and results of the Born cross section measurement for e*e^ — nhe. The 
table includes the integrated luminosity £, the number of observed signal events Nops, the radiative 
correction (1 + ô) and vacuum polarization correction factor |1 + II|?, the sum of the products of 
the branching fraction and efficiency Y? €;B;, the Born cross section o? 
90% C.L.), and the statistical significance S. 


vs (MeV) L (pb) Nas (14-9) [1 - SZ ef (1077) c? (pb) S 


and its upper limit (at the 


4085.4 524 0.07)? 0.68 1.052 2.40 0.0724 + 5.4 (< 23.7) 0.0c 
4188.6 43.1  0.0'2? 0.69 1.056 2.24 0.0120-6 + 13.7 (< 52.2) 0.00 
4207.7 546 251 0.75 1.057 2.22 21.81125 E57 (< 53.6) (ré 
4217.1 54.1 0.8729 0.85 1.057 2.18 3.8124 + 1.0 (< 32.2) 0.50 
4226.3 1091.7 41.273? 0.95 1.056 1.97 9.5120 E 2.7 5.80 
4241.7 55.6 0.07)? 1.06 1.056 1.72 0,07 59.3: 5.0 (< 17.6) 0.0c 
4258.0 825.7 10.3728 1.11 1.054 1.56 E + 1.2 (< 8.3) 2.00 
4307.9 449 0.0727 0.93 1.052 1.80 0.0*1 9 + 8.4 (< 35.3) 0.00 
4358.3 539.8 19.0729 0.81 1.051 2.07 10.0137 + 2.6 (< 19.3) 4.30 
4387.4 55.2 0.072? 0.90 1.051 1.87 0.070 + 5.8 (< 26.2) 0.00 
4415.6 1073.6 18.6775 0.94 1.053 1.65 5.3122 + 1.4 (< 11.2) 2.9 0 
4467.1 109.9 3.1'21 0.85 1.055 1.79 8.8183 + 2.3 (< 19.0) llo 
4527.1 110.0 2.1433 0.94 1.055 1.38 10*.2 + L8 (< 27.7) 0.80 
45745. 477 OOF? 115 1.055 0.88 0.0*115 + 6.8 (< 28.6) 0.00 
4599.5 566.9 4.0'23 1.27 1.055 0.75 35 2 Dg (< 11.1) 1.70 


summing them in quadrature. The following subsections describe the procedures and as- 
sumptions that led to these estimates of the uncertainties. 


A. Luminosity 


The integrated luminosity is measured using Bhabha events, with an uncertainty of 


1.0% [38]. 


B. Signal shape 


In the fit procedure, a discrepancy in the mass resolution between data and MC, as well 
as choices of background shapes and fit range introduce uncertainties on the results. Since 
the statistical fluctuation is large in the data sets, we cannot obtain a stable and reasonable 
estimation by simply comparing two fits with different choices. To avoid the influence of 
statistical fluctuations, ensembles of simulated data samples (toy MC samples) are generated 
according to an alternative fit model with the same statistics as data, then fitted by the 
nominal model and the alternative model. These trials are performed 500 times, and the 
deviation of mean values in the two trials is taken as the systematic uncertainty. The data 
samples taken at vs = 4.226, 4.258, 4.358, and 4.416 GeV are used to obtain an average 
uncertainty. 

A discrepancy in mass resolution and mass scale between data and MC simulation affects 
the fit result. To estimate this uncertainty, the signal shape is smeared and shifted by con- 
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volving it with a Gaussian function with a mean value of —1.2 MeV and standard deviation 
of 0.04 MeV, which are obtained from the study of a control sample of e*e^ > nJ/v. Toy 
MC samples are generated according to the smeared MC shape and fitted with a smeared 
and unsmeared signal shape. The average deviation determined from the four data samples 
is 7.596 and is taken as systematic uncertainty. 


C. Background shape 


Similarly, to estimate the uncertainty due to the background shape, a sum of signal shape 
and a second-order polynomial function with parameters determined from the fit on data is 
used to generate toy MC, then the toy MC samples are fitted by models with a first-order 
and a second-order polynomial background, respectively. The average deviation from the 
four data samples is found to be 6.3% and is taken as systematic uncertainty. 


D. Fitting range 


The systematic uncertainty for the fit range is determined by varying the fit ranges 
randomly for 400 times. The standard deviation of the fit results is taken as systematic 
uncertainty, which is determined to be 2.8% from the four data samples. 


E. B(he > ync)B(n > y7) 


The branching fraction of he + yn. is taken from Ref. [53]. The uncertainty in this 
measurement is 15.7% and the uncertainty of Bin — yy) is 0.5% 45]. These uncertainties 
propagate to the cross section measurement. 


F. ISR correction 


To obtain the ISR correction factor, the energy dependent cross section is parameterized 
with the sum of 3 coherent BW functions fitted to the cross sections measured in this analysis 
and the CLEO value at 4.17 GeV [a4]. The uncertainty of the input cross section is estimated 
by two alternative models. First, the energy-dependent cross sections are fitted with a sum 
of BW and a second order polynomial function. Second, the cross sections are fitted with 
a second order polynomial function only. The maximum difference in ISR correction factor 
and detection efficiency among these hypotheses is taken as systematic uncertainty due to 
the ISR correction. 


G. Vacuum polarization correction 
To investigate the uncertainty due to the vacuum polarization factor, we use two available 


VP parameterisations [51] [4]. The difference between them is 0.3% and is taken as the 
systematic uncertainty. 
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H. 5 ec Bin, =>} Xi) 


The branching ratios B(n, — X;) are taken from BESIII measurements |55], and the 
uncertainty of each channel is given in Table The systematic uncertainties associated 
with the efficiency include many items: tracking, photon and PID efficiency, K}, n°, o and 
Ne reconstruction, kinematic fit, cross feed and size of the MC sample. The procedure to 
estimate each item is described below, and the results are also listed in Table 


e Charged track, photon reconstruction and PID efficiencies 


Both the tracking and PID efficiency uncertainties for charged tracks from the in- 
teraction point are determined to be 1% per track, using the control samples of 
Ji > nn 095, J/ > ppr and J/ > KLK ccc. . The uncertainty due 
to the reconstruction of photons is 1% per photon and it is determined from studies 
of e*e- — yuu control samples [57]. 


e K? efficiency 
The uncertainty caused by K$ reconstruction is studied with the processes J/) > 
K**K* and Ji > óK$K-*m*. The discrepancy of KẸ reconstruction efficiency 
between data and MC simulation is found to be 1.296 and is taken as systematic 
uncertainty. 


e n/n’ efficiency 


To estimate the uncertainty due to the resolution difference in M (^y) between data 
and MC simulation in the 7 and 7° candidate selection, the MC shape of 7 (1?) is 
smeared by convolving it with a Gaussian function that represents the discrepancy 
of resolution and is determined by the study of an ete~ — nJ/w control sample. 
The difference of reconstruction efficiencies with and without smearing is taken as 
systematic uncertainty. 


e ne decay model 


We use phase space to simulate 7. decays in our analysis. To estimate the system- 
atic uncertainty due to neglecting intermediate states in these decays, we study the 
intermediate states in 7. decays from (3686) > zm. Ne > X; and generate MC sam- 
ples accordingly. For channels with well-understood intermediate states, MC samples 
with these intermediate states are generated according to the relative branching ratios 
given by PDG [45]. The spreads of the efficiencies obtained from the phase-space and 
alternative MC samples are taken as the systematic uncertainties. 


e n line shape 


The uncertainties of the 7. line shape originate from the model of ne and the errors of 
its resonant parameters. In the current MC generator, the 7. line shape is described 
by a BW function. However, in E1 transitions he — yn. a cubic photon energy term 
with a damping term at higher energies is introduced to the signal shape because of 
the transition matrix element and phase space factor. To estimate this uncertainty, 
toy MC samples, generated according to the model that takes the E1 photon energy 
dependency into account, are analyzed to obtain the efficiency difference. The uncer- 
tainties due to the 7, resonant parameters are considered by varying My. and I, in 
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the MC simulation within their errors given by PDG . The sum of these two items 
added in quadrature is taken as systematic uncertainty due to the 7, line shape. 


e Kinematic fit 


For the signal MC samples, corrections to the track helix parameters and the corre- 
sponding covariance matrix for all charged tracks are made to obtain improved agree- 
ment between data and MC simulation [58]. The difference between the obtained 
efficiencies with and without this correction is taken as the systematic uncertainty due 
to the kinematic fit. 


e Cross feed 


To check the contamination among the 16 decay modes of ne, 40,000 MC events for 
each channel are used to test the event misjudgement. 


e Size of the MC sample 


The efficiency of each channel is obtained by MC simulation. The statistical uncer- 
tainty is calculated according to a binomial distribution. 


In the fit procedure, ejB(n. > ASS ein, — X;) is used to constrain the strength 
among different ne decay modes, so the uncertainty from e;jB(y. — X;) will affect the fit 
results. In this case, we cannot simply add the uncertainty from e Bin, — X;) in quadrature 
with the other uncertainties. To consider the uncertainties of e;B(y. — X;) and their 
influence to the simultaneous fit, we change the e;B(ry — X;) within their errors and refit 
the data sample. The change of the cross section with the new results is taken as systematic 
uncertainty. 

In this procedure, systematic uncertainties are divided into two categories: the corre- 
lated part, which includes tracking, photon efficiency, PID efficiency, 7°/n/Kg efficiency, ne 
line shape and kinematic fit, and the uncorrelated part, which includes the 7, decay mode, 
cross feed, MC samples and Bi, + X). These uncertainties are assumed to be distributed 
according to a Gaussian distribution. The uncertainties of the correlated part are changed 
dependently (increasing or decreasing at the same time for all channels), while the uncer- 
tainties of the uncorrelated part are changed independently. We change the uncertainties 
(both correlated and uncorrelated parts) with a Gaussian constraint and refit the data set 
500 times. The cross sections calculated with these trials are fitted with a Gaussian func- 
tion, whose standard deviation is taken as systematic uncertainty. To obtain a conservartive 
estimation, the maximum deviation of 16.7% from the data samples at vs = 4.226, 4.258, 
4.358 and 4.416 GeV is adopted as systematic uncertainty from » 7; e; x Bi, — X;) for all 
the data sets. 


VII. UPPER LIMIT WITH SYSTEMATIC UNCERTAINTY 


For the data sets without significant nh, signals observed, an upper limit at the 90% 
C.L. on the cross section is set using a Bayesian method, assuming a flat prior in c. In 
this method, the probability density function of the measured cross section c, P(c), is 
determined using a maximum likelihood fit. The 9096 confidence limit (L) is then calculated 
by solving the equation 
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TABLE II: Summary of systematic uncertainties on cP(e*e- — nh.) (in 96) at /s = 4.226 GeV. 


Sources uncertainty in e 
A. Luminosity 1.0 
B. Signal shape 7.5 
C. Background shape 6.3 
D. Fitting range 2.8 
E. Blhe > *n)B(n > y7) 15.7 
F. ISR correction 13.9 
G. VP correction 0.3 
H. Miei x Bine => X4) 16.7 
Total 28.7 


TABLE III: Systematic uncertainties (in 96) for e Bin, + X;) for each ne exclusive decay channel. 


S mune pp "n 2(K^K-) isa Dpnm- "us K*'K 2(n'm-) nme 
Tracking eff. 2.0 4.0 4.0 4.0 4.0 6.0 6.0 2.0 
Photon eff. 3.0 3.0 3.0 3.0 3.0 3.0 3.0 5.0 
PID 2.0 4.0 4.0 4.0 4.0 6.0 6.0 2.0 
K? eff. — = — ES - — — — 

n eff. x = = — - — — 3.0 

n eff. 1.7 2.0 2.6 1.7 p 1.7 2.0 1.3 

Ne decay model 0.0 2.1 3.7 0.6 2.5 0.0 3.0 4.6 

ne line shape 6.2 6.2 6.2 6.2 6.2 6.2 5.0 5.0 
Kinematic fit 2.3 3.8 3.9 3.5 3.0 6.1 4.4 1.3 
Cross feed 0.0 0.7 0.0 2.0 0.0 0.0 0.0 0.0 
MC sample 0.9 1.2 1.0 1.4 1.4 1.6 2.1 1.4 
B(n.— Xi) 37.0 220 46.0 26.0 34.0 — 280 54.0 23.0 

oo 
0.1 = f P(o)do. (3) 
L 


To include multiplicative systematics, Pio) is convolved with a probability distribution 
function of sensitivity, which refers to the denominator of Eq. and is assumed to be a 
Gaussian with central value S and standard deviation o, [59]: 


Pe J "P (Šo) T E ds. (4) 


Here, P(c) is the likelihood distribution obtained from the fit and parameterized as double 
Gaussian. By integrating P'(c) we obtain the 90% C.L. upper limit taking the systematic 
uncertainties into account. 
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Sources PPT ee EE SEET DEE EARN EI ET E EE EI EEN 
Tracking eff. 2.0 4.0 6.0 2.0 2.0 4.0 2.0 4.0 
Photon eff. 5.0 3.0 3.0 5.0 5.0 5.0 7.0 7.0 
PID 2.0 4.0 6.0 2.0 2.0 4.0 2.0 4.0 
K? eff. = 12 1.2 - - - - - 
T? eff. 2.3 - - - - - 2d 1.5 
1) eff. 2.2 1.6 2.1L 1.7 1.5 2.2 1.9 1.2 
Ne decay model 5.8 2.5 5.2 5.5 8.1 0.0 0.1 0.5 
Nc line shape 5.1 5.0 6.2 5.0 6.2 5.0 5.1 5.1 
Kinematic fit 2.1 2.4 L7 1.1 2.5 1.8 0.5 2.9 
Cross feed 0.0 0.0 0.0 0.0 7.7 0.0 0.0 0.0 
MC sample 1.4 1.4 2:1 1.4 1.5 1.9 1.9 2.9 
B(n.— Xj) 380 210 28.0 280 540 ` 300 22.0 20.0 


VIII. RESULTS AND DISCUSSION 


In this study, the Born cross section and its upper limits of e^ e^ — nh, are measured with 
statistical and systematical uncertainties at c.m. energies from 4.085 to 4.600 GeV, and the 
results are listed in Table [| Clear signals of e*e^ — nhe are observed at vs = 4.226 GeV for 
the first time. The Born cross section is measured to be (9.5732 + 2.7) pb. We also observe 
evidence for the signal process at va = 4.358 GeV with a cross section of (10.0737 + 2.6) pb. 
For the other c.m. energies considered, no significant signals are found, and upper limits 
on the cross section at the 90% C.L. are determined. The cross sections measured in this 
analysis and CLEO [54 are modeled with a coherent sum of three BW functions (as shown 
in Fig. [5) to calculate the ISR correction factors. 

Comparing with the process ete~ > nJ/v 60 , if we suppose both processes come from 
higher mass vector charmonia, the ratio r(Y — nh.)/T(W — nJ/w) is determined to be 
0.20 + 0.07 and 1.79 + 0.84 at 4/s = 4.23 GeV and 4.36 GeV, respectively. These results are 
larger than theoretical expectation: l'((4160) > gh) /T(v(4160) > nJ/w) = 0.07887 and 
Dl(v(4415) > qh.)/T(v(4415) > nJ/) = 0.06736 

Comparing with the cross section of ee > ntn he [33], we find that the cross section 
of ee — nh, is smaller. But due to the limited statistics we cannot determine the line 
shape of c.m. energy dependent cross section precisely. 
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